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Moment analysis of the probability distributions of different sandpile models 
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We reconsider the moment analysis of the Bak-Tang-Wiesenfeld and the Manna sandpile model 
in two and three dimensions. In contrast to recently performed investigations our analysis turns 
out that the models are characterized by different scaling behavior, i.e., they belong to different 
universality classes. 



I. INTRODUCTION 

The Bak-Tang-Wiesenfeld (BTW) model was intro- 
duced as a paradigm of the concept of self-organized 
criticality which describes the emergence of spatiotempo- 
ral correlations in slowly driven dissipative systems [0,0 . 
Despite its analytical tractability |^ the scaling behav- 
ior of the two-dimensional BTW model is not well un- 
derstood. Especially the exponents which determines 
the avalanche distributions are not known exactly. Sev- 
eral numerical attempts were made but do not provide 
consisting results p|-[l()[ . Recently De Menech et al. per- 
formed a moment analysis of the BTW model |ll| which 
was extended by several authors to different sandpile 
models fl^-Q. Especially the moment analysis of the 
size distribution of the BTW and Manna sandpile model 
has led Chessa et al. to the conclusion that both models 
are characterized by the same scaling exponents and thus 
belong to the same universality class |12|. In this work 
we reconsider the moment analysis and compare the scal- 
ing behavior of various avalanche quantities for the BTW 
and Manna model. Our analysis turns out that in con- 
trast to [Q the moment behavior of both models differs 
significantly, i.e., the BTW and the Manna model belong 
to different universality classes. 



II. MODELS AND SIMULATIONS 

The BTW model is defined on a D-dimensional square 
lattice of linear size L in which non-negative integer vari- 
ables Ej. represent a physical quantity such as the lo- 
cal energy, stress, height of a sand column, etc. One 
perturbes the system by adding particles at a randomly 
chosen site r according to 



E. 



Er. 



1. 



with random r. 



(1) 



A site is called unstable if the corresponding variable ex- 
ceeds a critical value Ec, i.e., if E^ > Ec, where the 
critical value is given by Ec — 2D. An unstable site re- 
laxes, its value is decreased by Ec and the 2D nearest 
neighboring sites are increased by one unit, i.e.. 



E, 



Er, 



Er- 



En 



Ec 



+ 1. 



(2) 



(3) 



In this way the neighboring sites may be activated and 
an avalanche of relaxation events may take place. These 
avalanches are characterized by several physical proper- 
ties like the size s (number of relaxation events) , the area 
a (number of distinct toppled sites), the time t (num- 
ber of parallel updates until the configuration is stable), 
the radius r (radius of gyration), the perimeter p (num- 
ber of boundary sites), etc. In the critical steady state 
the corresponding probability distributions should obey 
power-law behavior [ll 



Px{x) 



(4) 



characterized by the avalanche exponents Tx with x G 
{s^ajt^r^p}. Assuming that the size, area, etc. scale as 
power of each other, 

X ~ x'''""' (5) 

one obtains the scaling relations ^xx' — (jx' — \)/ {jx — 1). 
The scaling exponents ^xx' describe the static avalanche 
properties as well as its propagation. For instance, the 
exponent ^sa indicates if multiple toppling events are rel- 
evant (7sa > 1) or irrelevant {"fsa = 1)- The exponent 
7ar equals the fractal dimension of the avalanches. A 
possible fractal behavior of the avalanche boundary cor- 
responds to the inequality D ~ \ < jpr < D. Finally, 
the exponent jtr is usually identified with the dynamical 
exponent z. 

A stochastic version of the BTW model was introduced 
by Manna ||l^. Here, critical sites relax to zero, i.e., 
Er — > if Er > Ec and the removed energy is randomly 
distributed to the nearest neighbors in the way that one 
chooses randomly for each energy unit one neighbor. For 
Ec = 1 the behavior of the model corresponds to a simple 
random walk. Above this value {Ec > 2) is the choice of 
the critical energy irrelevant to the scaling behavior (see 

Fig.0). 

Recently Dhar introduced a modified version of the 
two-dimensional Manna model where the energy of crit- 
ical sites is not reduced to zero but Eij -^ Eij — 2. The 
energy AE = 2 is then equally distributed with proba- 
bility 1/2 to the sites {i ± 1, j) or otherwise to the sites 
(«, j±l) ||l^. In this case it is possible to extend an oper- 
ator algebra, which was successfully applied in studying 
the BTW model ||^, to this modified Manna model. 

Compared to the BTW model the dynamics of the 
Manna model with its stochastic distribution of the en- 
ergy to the nearest neighbors can be interpreted as a 
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FIG. 1. The probability distribution Pa{a) of the 
two-dimensional Manna model for E^ G {2, 3, 5, 10} and vari- 
ous system sizes L. For L < 512 the curves are shifted in the 
downward direction. No significant difference between the 
curves for different values of the critical energy can be ob- 
served, i.e., Ec does not affect the scaling behavior of model. 



disorder effect. A different kind of disorder effects were 
investigated in directed sandpile models by introducing 
stochastic toppling conditions. In particular the exact 
solution of the directed BTW model [0| was used as a 
starting-point in order to examine how the known scaling 
behavior of the system is affected by the stochastic top- 
pling rules. Depending on the details of the introduced 
disorder a collapse of the critical behavior fl^, non- 
universal critical behavior pSpG] as well as a crossover to 
the different universality class of directed percolation ||2l| ] 
was observed. 

Thus an important question is if the scaling behavior 
of the Manna model differs from the scaling behavior of 
the BTW model, i.e., if the additional fluctuations of the 
Manna model could change the universality class. A real- 
space renormalization scheme predicted that both models 
belong to the same universality class [p2|,p3|. Here, the 
authors used a mean-field-type approximation p4| in or- 
der to perform a block transformation. Therefore it is 
not clear if this renormalization ansatz is an appropri- 
ate tool to take the additional fluctuations of the energy 
distribution of the Manna model into account. 

A momentum space analysis of Langevin equations of 
the BTW and the related Zhang model pS] predicted 
that both models are characterized by the same expo- 
nents p6lE7| and numerical investigations confirmed this 
prediction|l8|,p8| . Unfortunately up to now no extension 
of the momentum space analysis to the Manna model 
could be performed. The crucial point of this renormal- 
ization approach is the choice of an appropriate noise- 
correlator 12^. Compared to the BTW model the corre- 
sponding analysis of the Manna model requires a different 
noise-correlator. But a different noise-correlators could 
lead to a different critical behavior (see for instance P9|). 

The first numerical indication that both models be- 
long to different universality classes was reported by Ben- 



Hur and Biham g . They measured several scaling expo- 
nents jxx' via Eq. (^) and found, in particular, that the 
dynamical exponents z, the multi-toppling exponents 7sa 
and the boundary exponents jpr of both models differ 
significantly. But one has to mention that at least one 
of the scaling exponents of the BTW model (7sa) can- 
not be determined in this way since the assumed scaling 
behavior s ~ a'''"" is not well defined |^ 12|. 

However, the conjecture of Ben-Hur and Biham was 
confirmed for D = 2 by a numerical determination of the 
avalanche exponents Ts, Ta, Tt, and Tr W which again 
differ significantly for both models. In this analysis the 
exponents of the Manna model were obtained by a direct 
regression analysis as well as a simple finite-size scaling 
analysis of the corresponding probability distributions. 
In the case of the BTW model it was found that the 
probability distributions are affected by unconventional 
logarithmic finite-size corrections MM which lead to un- 
certain results for the simple finite-size scaling ansatz. 
Taking these corrections into account it is possible to 
estimate the values of the exponents |Q,y by an extrapo- 
lation to L ^ oo. But one has to note that the assumed 
logarithmic corrections are found only numerically, i.e., 
there exist up to now no analytical justification of this 
unconventional behavior. 

These difficulties vanish in the three-dimensional case 
where both models fulfill the finite-size scaling ansatz 
which makes the analysis much easier JS^]. The accu- 
racy of the determination is sufficient to show that both 
models belong to different universality classes [Q. Ad- 
ditionally the scaling behavior of the three-dimensional 
Manna model is strongly affected by multiple toppling 
events (ts ^ Ta) whereas it seems that the rare mul- 
tiple toppling events of the BTW model does not con- 
tribute to the scaling behavior |||,^ . Taking this results 
into consideration we have convincing, but of course not 
completely rigorous, arguments that the BTW and the 
Manna model does not belong to the same universality 
class. 

Recently this statement was questioned by Chessa et 
al. ||l^ who performed a moment analysis of both models 
in analogy to the investigations of De Menech et al. [11| . 
The g-moment of the probability distribution Px{x) is 
defined as 



{x'i 



dx x'^ Px{x). 



(6) 



The finite system size L causes a cut-off of the probabil- 
ity distribution at Xmax ~ L'^'''^ . Assuming a power-law 
behavior of the distributions [Eq. (0)] the scaling behav- 
ior of the (/-moment is dominated by the upper boundary 
of Eq. (i) 



{^')l 



Aq) 



(7) 



if g > Ta; — 1 and where the moment exponent Uxici) is 
given by 



(^x{q) = lxr{q+ 1 -Trr). 



(8) 



For q < Tx ~ 1 the moment exponent (Jx{<l) behaves 
nonlinear with respect to q. The normalization of the 
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FIG. 2. The exponent (^aisi) and 

(7r(g) of the two-dimensional Manna model. The extrapo- 
lation to the horizontal axis yields the exponent Ts = 1.286 
and Tr — 1.729, respectively. The slope corresponds to the 
scaling exponents 7sr and 7rr where the latter equals one 
(see Table M). The positions and sizes of the circles on the 
horizontal axes correspond to the values and error-bars of the 
exponents Ts = 1.275±0.011 and r^ = 1.743±0.025, obtained 
from a regression analysis M. 

probability distributions results in Uxiq) — for g — > 0. 
Performing numerical investigations one can obtain the 
behavior of the exponent ax (q) via a regression analysis 
of ln(x^}L as a function of Ini, or via the logarithmic 
derivative 



(^x{q) 



dlnL 



(9) 



The behavior of ax (q) as a function of q allows to de- 
termine the scaling exponent "fxr, which corresponds to 
the slope, and the avalanche exponent Tx which can be 
obtained from an extrapolation to the horizontal axis 
[Eq. (^)]. This is shown in Fig. g where the expo- 
nents as (q) and ar (q) of the Manna model are plotted as 
a function of q. The values of the exponents are in agree- 
ment with those of previous investigations obtained from 
a regression analysis and a finite-size scaling analysis, re- 
spectively [H . But one has to mention that the accuracy 
of the regression analysis is higher than the accuracy of 
the moment analysis. In the case of the regression anal- 
ysis one obtains the exponents Tx by a direct fit to the 
distribution. Whereas for the moment analysis one has 
first to calculate the average [Eq. (pi)], second to fit the 
logarithmic derivatives [Eq. (||)] and third to extrapolate 
to ax{q) = 0. This leads to a propagation of errors which 
increases the uncertainty significantly. 

The g-dependence of the moment exponent ax{q) is 
determined by the avalanche exponents [Eq. (||)]. There- 
fore, the moment analysis can be used to distinguish the 
universality classes of different models. Using the scaling 
relation '^xr ~ {Tr — i)/{Tx — 1), Eq. (j^ now reads 



f^xiq) ^Ixrq 



(10) 



with S = 1 — r,.. Thus, we get that the intercept E of 



the linear ^-dependence of the moment exponent ax (q) is 
the same for all distributions (site, area, duration, etc.) 
and is therefore a characteristic quantity of the model. 
Considering two models we get that different values of 
S implies different universality classes. But the same 
value of the intercept S does not imply that both models 
belong to the same universality class. It is possible that 
the models display different values of 'jxr which results 
in different values of r^: = 1 — "E/jxr- Thus two different 
models belong to the same universality class if they are 
characterized by the same linear dependence of ax (q) for 
all relevant quantities x. 

The determination of the intercept E for various dis- 
tributions allows to estimate the accuracy of the mo- 
ment analysis. In the case of the Manna model we ob- 
tain from the moment analysis of the size, area, duration 
and radius distribution the values E^ — —0.7900 ±0.002, 
Ea = -0.7202 ± 0.003, Ej = -0.7684 ± 0.004, E^ = 
—0.7333 ± 0.005. The corresponding values of the scal- 
ing exponents ^xr are listed in Table |. The above 
error-bars correspond to the uncertainty of the linear 
regression [Eq. (iffl]. The average value EManna.2D — 
0.7530 ± 0.037 agrees with E = 3/4 predicted in [§. 
The latter error-bar reflects the uncertainty of the whole 
method and could be use as a lower bound for the error 
of the avalanche exponent Ar^ > 0.037. Typical error- 
bars of a direct analysis of the probability distributions 
are Ar^ < 0.025. 



III. MOMENT ANALYSIS OF THE BTW AND 
MANNA MODEL 

We now consider the moment analysis of the BTW 
and Manna model for the size, area, duration, radius, 
and perimeter distribution and compare the correspond- 
ing results. In the case of the BTW model we analyze the 
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FIG. 3. The exponent a3{q) of the two-dimensional BTW 
and Manna model. The curves differs slightly. The derivative 
of both curves is plotted in the inset and reveals that the 
exponent as (q) of the BTW model is not well defined in the 
sense that no clear saturation of the exponent can be seen for 
q>l. 



accepted for publication in Physical Review E. 



probability distributions of lattice sizes up to L = 4096 
ior D = 2 and L — 256 for D = 3, respectively. Since the 
Manna model does not display strong finite-size effects as 
the BTW model it is sufficient to consider system sizes 
only up to i == 2048 for Z? = 2. 

The moment analysis of the size distribution leaded 
Chessa et al. to the conclusion that both models are char- 
acterized by the same scaling behavior [ 121 . Therefore we 
first consider the size distribution and the correspond- 
ing exponent as{q) is shown in Fig. g[ In contrast to 
Chessa et al. who found that both models display in- 
distinguishable curves for g > 1 we get slightly different 
curves. Plotting the derivative of the exponents as{q) 
with respect to q the differences become more signifi- 
cant (see inset of Fig. 0). In the case of the Manna model 
the derivative of the exponent saturates quickly with in- 
creasing q. Whereas the derivative dqas{q) for the BTW 
model is characterized by a finite curvature, i.e., the ex- 
ponent 'Ysr is not well defined in the case of the BTW 
model. The duration probability distribution displays a 
similar behavior as it can be seen in Fig. H. Again the 
exponent at{q) of the BTW model is characterized by a 
finite curvature. 

The cr-exponents of the radius and area distribution 
are plotted for both models in Fig. H. In the case of the 
BTW model the area and radius distributions display in 
contrast to the size and duration distribution the usual 
behavior. The g-dependence of the exponents Ua{q) and 
(Tr(a) is given by a linear function and the slopes corre- 
spond to the trivial values ^ar — 2 (compact avalanches) 
and 7rr = 1 (see Table |) . These exponents are the same 
for the Manna model and the corresponding curves are 
parallel. But as Fig. g shows there is a clear shift be- 
tween the curves of the BTW and Manna model, i.e., 
both models are characterized by two different avalanche 
exponents Tq and r^. 

The analysis of the perimeter distribution yields again 
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FIG. 4. The exponent <Jt{q) of the two-dimensional BTW 
and Manna model. The derivative of both curves is plotted 
in the inset and reveals that the exponent at (g) of the BTW 
model displays the same complicated behavior as the expo- 
nent (Js{q). 
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FIG. 5. The exponent Caiq) and (Tr(g) (inset) of the 
two-dimensional BTW and Manna model. One obtains for 
both models the same asymptotic slopes (see Table |y) which 
agree with 7ar = 2 and -yrr ~ 1. But the extrapolation to the 
horizontal axis yields different values for the models indicat- 
ing that the BTW and the Manna model are characterized by 
different avalanche exponents Ta and r^, respectively. 



a simple linear g-dependence of the exponents ap (q) (not 
shown) and we obtain 7pr = 1.266 ± 0.019 which corre- 
sponds to the fractal dimension of the boundary. This 
value is in agreement with that obtained from a direct 
analysis of the scaling relation p ~ r'^'"' and differs signif- 
icantly from the value of the Manna model 7pr = 1.42 [B|. 

Due to the nonlinear behavior of the size and dura- 
tion distribution we use the area, radius, and perime- 
ter distributions in order to estimate the intercept E of 
the BTW model and obtain Ebtw.2D = 0.391 ± 0.011. 
But one has to be carefully to compare this result with 
the corresponding value of the Manna model. The per- 
formed moment analysis based on the assumption that 
the scaling behavior of the BTW model is given by a 
pure power-law behavior [Eq. (0)]. But this assumption 
is in contradiction to the observed logarithmic finite-size 
corrections MM as well as to recently reported investi- 
gations where a multifractal behavior of the distribution 
was observed O]. In the latter case no avalanche expo- 
nent could be defined whereas the analysis of the loga- 
rithmic finite-size corrections yield Sbtw,2£i ~ 2/3 M 
which again differs significantly from the corresponding 
value of the Manna model. 

These results lead us to the conclusion that the mo- 
ment analysis is applicable in order to determine the ge- 
ometric properties of the avalanches of the BTW model, 
e.g. 7ar, Jpr- The obtained results are in agreement with 
those of previous investigations which based on different 
methods of analyzing. But in contrast to the geometric 
properties the moment analysis of the dynamical proper- 
ties (size and duration) of the avalanches exhibit a non- 
trivial behavior. This shows that the assumed simple 
power-law behavior is not fulfilled and that scaling cor- 
rections to Eq. (^ cannot be neglected in these cases. 
The question whether these corrections can be inter- 
preted in terms of unconventional finite-size effects [^,^1 , 
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FIG. 6. The exponent f7a{q) of the three-dimensional BTW 
and Manna model. One obtains for both models the same 
asymptotic slopes (see Table |) which agree with ■jar = 3 
(compact avalanches). In the case of the BTW model both 
exponents crs{q) and Ua{q) are plotted, but no difference be- 
tween the curves can be seen. The inset shows the exponents 
o's{q) (solid line) and (Ja{q) (dotted line) for the Manna model. 
The different slopes indicate that multiple toppling events are 
relevant in the case of the Manna model {jsa = 73r/7ar > !)• 



or correspond to a cross-over effect from the boundary to 
the bulk regime R,|l^l , or indicate a multifractal behavior 
of the two-dimensional model [0| , remains open. 

In the following we briefly compare the three- 
dimensional BTW and Manna model and focus our at- 
tention on the analysis of the size and area probability 
distribution. In Fig. O we plot the moment exponent 
(Ja{<l) for both models. The result is similar to the two- 
dimensional case (Fig. a). The slopes of the curves agree 
with the value "far = 3 indicating that both models are 
characterized by compact avalanches (see Table |). But 
the two curves are shifted which shows that the avalanche 
exponents r^ of the BTW {Ta = 1.352±0.022) and Manna 
{Ta = 1.436 ± 0.018) model are different for D = 3. 

The determination of the intercepts E confirms this 
result. Analyzing for both models the size, area, duration 
and radius distribution we obtain Sbtw,3D — 1-016 ± 
0.056 and SManna,3D = 1-333 ± 0.036- The intercepts 
differs significantly. 



TABLE I. Exponents 7^,^ for the BTW and Manna model. 
The values are obtained from a linear regression according 
to Eq. (nol) and the errors are of the order Ajxr < 0.01. 
The real error which includes the uncertainty of the whole 
moment analysis is hard to estimate and increases this value 
significantly (see text). 





BTW, 2D 


Manna, 2D 


BTW, 3D 


Manna, 3D 


7sr 




2.764 


3.004 


3.302 


7ar 


2.021 


2.025 


3.004 


3.076 


7tr 




1.540 


1.618 


1.713 


7pr 


1.266 


1.42 § 






7rr 


1.008 


1.006 


0.995 


1.010 



Finally we consider the area and size distribution for 
both models. As mentioned above the three-dimensional 
Manna model behaves different than the BTW model in 
the sense that multiple toppling events affects the scaling 
behavior decisively. This result can also be obtained from 
the moment analysis. In the case of the BTW model no 
significant difference between the exponents Uaiq) and 
(Js{q) can be observed (see Fig. g). Whereas for the 
Manna model both curves are characterized by different 
slopes as can be seen in the inset of Fig. g. A regression 
analysis yield 7^^ = 3.302 ± 0.06 and "far = 3.076 ± 0.09, 
respectively. Thus multiple toppling events affect the 
scaling behavior of the Manna model strongly (r^ 7^ r^), 
in contrast to the BTW model. 



IV. SUMMARY 

In summary we performed a moment analysis of several 
probability distributions of the BTW and Manna model 
for D = 2 and D = 3. We found that the correspond- 
ing moment exponents (Jx{q) differ significantly for both 
models showing that the BTW and Manna model belong 
to different universality classes. Recently performed sim- 
ulations of sandpile models on a Sierpinski gasket dis- 
pla y agai n the different scaling behavior of both mod- 
els 



31 



33]. Our results confirm the universality hypoth- 
esis of Ben-Hur and Biham where the scaling behavior 
of sandpiles models is determined by the way in which 
the relaxing energy of critical sites is distributed to the 
neighboring sites ||g]. 
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